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If dark matter is embedded in a non-trivial dark sector, it may annihilate and decay to lighter 
dark-sector states which subsequently decay to the Standard Model. Such scenarios - with annihi¬ 
lation followed by cascading dark-sector decays - can explain the apparent excess GeV gamma-rays 
identified in the central Milky Way, while evading bounds from dark matter direct detection experi¬ 
ments. Each ‘step’ in the cascade will modify the observable signatures of dark matter annihilation 
and decay, shifting the resulting photons and other final state particles to lower energies and broad¬ 
ening their spectra. We explore, in a model-independent way, the effect of multi-step dark-sector 
cascades on the preferred regions of parameter space to explain the GeV excess. We find that the 
broadening effects of multi-step cascades can admit final states dominated by particles that would 
usually produce too sharply peaked photon spectra; in general, if the cascades are hierarchical (each 
particle decays to substantially lighter particles), the preferred mass range for the dark matter is 
in all cases 20-150 GeV. Decay chains that have nearly-degenerate steps, where the products are 
close to half the mass of the progenitor, can admit much higher DM masses. We map out the region 
of mass/cross-section parameter space where cascades (degenerate, hierarchical or a combination) 
can fit the signal, for a range of final states. In the current work, we study multi-step cascades in 
the context of explaining the GeV excess, but many aspects of our results are general and can be 
extended to other applications. 

PACS numbers: 95.35.+d, 12.60.-i; MIT-CTP/4647 


I. INTRODUCTION 

Over the past five years, numerous independent stud¬ 
ies have confirmed a flux of few-GeV gamma rays from 
the inner Milky Way, steeply peaked toward the Galactic 
Center, that is not captured by models for the known 
diffuse backgrounds mm- This “Galactic Center ex¬ 
cess” (GCE), detected using public data from the Fermi 
Gamma-Ray Space Telescope, has a spatial morphology 
well described by the square of a generalized Navarro- 
Frenk-White (NFW) profile, projected along the line of 
sight. Furthermore, it is highly spherically symmetric, 
centered on the Galactic Center (GC), and extends at 
least 10 degrees from the GC [lOjQ these conclusions 
remain unchanged when accounting for systematic un¬ 
certainties in the modeling of the diffuse backgrounds 
m These spatial properties suggest the excess emission 
could arise from the annihilation of dark matter (DM) 
with an NFW-like density profile. Competing interpre¬ 
tations include a transient event at the GC producing 
high-energy cosmic rays that subsequently yield few-GeV 
gamma rays by scattering processes dun, or a popula¬ 
tion of many unresolved millisecond pulsars (MSPs) (e.g. 
mm)- However, these interpretations face significant 
challenges: it is unclear whether the proposed outflow 
models can match the spectrum and morphology of the 
excess [T(5] (see also mm), and estimates of the MSP 
population in the region of interest consistently under¬ 
predict the signal by an order of magnitude [19, 20). 


1 This analysis exploited improvements to the Fermi point spread 
function as described in [12- 


Models where DM annihilates with a roughly ther¬ 
mal cross-section and has a mass of order several tens 
of GeV can readily account for the spectrum and size of 
the excess. However, when embedded in even a simplified 
DM model, there are often powerful constraints on these 
scenarios from direct detection and collider bounds (e.g. 
eh na). While UV-complete models where the DM an¬ 
nihilates directly to Standard Model (SM) particles do 
exist (e.g. (25H25] ). the constraints are much more eas¬ 
ily evaded if the DM produces gamma-rays via a cascade 
process pfiHHT)] . In such scenarios, the DM is secluded in 
its own hidden dark sector, and first annihilates to other 
dark sector particles; these mediators subsequently decay 
into SM particles that produce gamma-raysj^] 

The presence of an intermediate step between DM an¬ 
nihilation and the production of SM particles broadens 
the spectrum of SM particles produced, and consequently 
also broadens the resulting gamma-ray spectrum, unless 
the mediator is degenerate in mass with either the DM or 
the total mass of the SM decay products. The gamma- 
ray multiplicity is increased by a factor of two, if each 
mediator decays into two SM particles, and the typical 
energy of the gamma-rays is reduced accordingly. Thus 
cascade models for the excess generically tend to accom¬ 
modate: 

• Higher DM masses, 

• Decays of the mediator to SM final states whose 
decays produce a more sharply peaked gamma-ray 
spectrum than favored by direct annihilation. 


2 Annihilation into the dark sector can also lead to a novel spa¬ 
tial distribution for the signal eh, but the GCE favors a cuspy 
morphology, so in this work we assume all decays are prompt. 
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In general, there may be more than one decay step 
within the dark sector; the dominant annihilation of the 
DM need not be to the lightest dark sector particle (e.g. 
[321133]). If couplings within the dark sector are stronger 
than couplings between the sectors, dark sector parti¬ 
cles will preferentially decay within the dark sector, with 
decays to the SM only occurring when no other states 
are available. Regardless of the model under consider¬ 
ation, in the absence of a mass degeneracy, each decay 
will increase the final gamma-ray multiplicity, decrease 
the typical gamma-ray energy, and broaden the spectrum 
(in the presence of a mass degeneracy only the first two 
effects will occur). Accordingly, long decay chains could 
potentially permit much heavier DM to explain the GCE, 
or favor decays to different SM states. In a sense, this 
description also characterizes the known decays of SM 
particles; final states whose decays produce gamma-rays 
through a lengthy cascade will generate a broader spec¬ 
trum with a lower-energy peak, compared to final states 
that generate gam ma-ra ys via a short cascade (we discuss 
this further in Sec. Ill). 


It is this possibility of multi-step dark sector cascades 
that we explore in this work. For simplicity, we con¬ 
sider the case where all dark-sector particles involved in 
the cascade (except possibly the DM itself) are scalars 
- we briefly discuss the case of non-scalar mediators in 
Sec. [TV] In this case, the results are largely independent 
of the details of the dark sector. The DM pair-annihilates 
into two scalar mediators which subsequently undergo a 
multi-step cascade in the dark sector, eventually produc¬ 
ing a dark-sector state (with high multiplicity) that de¬ 
cays to the SM: 


XX (frn'fin —> 2 X (j) n -i(f> n -i 


-m— 1 


x -> 2 n X //. 


( 1 ) 


Here // are SM lepton or quark pairs, which can subse¬ 
quently decay; the decays shown above may also produce 
photons in the final step via final state radiation (FSR). 
By fitting the resulting photon spectrum to the GCE, 
we determine the allowed values of cross-section and DM 
mass for cascades with one to six steps, for a variety of 
SM final states. Provided that the masses of the particles 
at each step in the cascade are not near-degenerate, the 
final spectrum of gamma-rays becomes nearly indepen¬ 
dent of the exact masses at each step. This assumption is 
not limiting, as results for the quality of fit for the more 
general case of non-hierarchical cascades (with nearly- 
degenerate steps) can be simply extracted from results 
derived assuming a large hierarchy. 


In Sec. |TT] we outline the determination of the photon 
spectrum for an n-step cascade with specified SM final 
state, and discuss the procedure used to compare such a 
spectrum to the GCE. We present sample results of these 


fits in Sec. Ill under certain assumptions. Section IV 


extends our results for general cascades, and contains 
our complete fit results. In Sec. [V] we outline the existing 
experimental constraints a complete model for the GCE 


via cascade decays would need to satisfy. We present 
our conclusions in Sec. |VI[ In the appendices we provide 
additional details of our methodology and discuss some 
further model-dependent considerations. 


II. METHODOLOGY 


The photon flux generated by the annihilations of self¬ 
conjugate dvQ as a function of the direction observed in 
the sky, is given by: 
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(av) dN 7 
87t dEj 


J {l,b) , 


( 2 ) 


where (av) is the thermally averaged annihilation cross- 
section, m x is the DM mass, and dN 1 /dE 1 is the photon 
spectrum per DM annihilation, which has contributions 
from FSR and from the decay of the leptons or quarks 
and their subsequent hadronization products. The J- 
factor, the integral of DM density squared along the line- 
of-sight, is a function of the observed direction in the sky 
expressed in terms of Galactic coordinates l and 6 : 


J(l,b) 




2 r Q scosZcos 6 



(3) 


where r Q « 8.5 kpc is the distance from the Sun to the 
Galactic Center, and s parametrizes the integral along 
the line-of-sight. We parameterize the DM density by a 
generalized NFW halo profile 04], 35]: 


P (r, 7 ) = Po 


(r/r s ) 7 
(1 + r/r s f “ 7 


(4) 


Here we use r s = 20 kpc, po = 0.4 GeV/cm 3 and 7 = 1.2, 
following m, as we will compare our models to the data 
using the spectrum and covariance matrix determined by 
that work. 

We focus on n-step cascades ending in cj) 1 —► //, where 
// is a pair of electrons, muons, taus or 6 -quarks. Other 
SM final states are possible, of course, but these cases 
span the range from steeply peaked photon spectra close 
to the DM mass through to the lower-energy and broader 
spectra characteristic of annihilation to hadrons. In order 
to generate the cascade spectrum, we first start with the 
result from direct DM annihilation, which is equivalent to 
the spectrum from (f>i decay (in the <f> 1 rest frame) if the 
DM mass is half the <j> 1 mass. For the case of electrons or 
muons we determine this spectrum analytically using the 
results of [36], whilst for taus and 6 -quarks the results are 
simulated in Pythia8 m- We have relegated the details 
of calculating these spectra to Appendix [A] 


3 As discussed in Appendix^] our results can be readily translated 
to the case of decays, although the steeply peaked morphology 
of the GCE disfavors this interpretation. 
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FIG. 1. Oth step (direct annihilation) photon spectra dN^/dxo for 
(j >i decaying to (e, /x, r, 6) in (blue, red, green, orange). Solid curves 
correspond to ef = 0.1, and dashed to ef = 0.3. The electron and 
muon spectra have been magnified by a factor of ten to appear 
comparable to the taus and 6s. 


We denote the spectrum obtained at this “Oth step” 
by dNry/dx o, where xq = 2Eo/mi , mi is the mass of 
0 i and Eq is the energy of the photon in the 0i rest 
frame. The shape of the photon spectrum is determined 
by the identity of the final state particle / and also the 
ratio £/ = 2mf Im\. In the limit where the decay of (f>i 
is dominated by a two-body final state (at least for the 
purposes of photon production), the photon spectrum 
converges to a constant shape (as a function of xq) as 
e/ —► 0 and the // become highly relativistic. However, 
final state radiation (FSR) and hadronization depend on 
the energy of the // products of the 0i decay in the 
rest frame, so in cases where these effects dominate, 
the dependence of the photon spectrum on ef is more 
complex. 

In Fig. [I] we show dN x /dx o per annihilation for the 
four different final states we considered, for e/ = 0.1 and 
ef = 0.3. The photon spectra from electron and muon 
production are dominated by FSR, whereas for 6-quarks 
fragmentation and hadronization are important. In the 
photon spectrum from taus, these effects are subdomi¬ 
nant and so the impact of varying e/ is minimal. Note 
that the spectrum for 6-quarks is peaked at a significantly 
lower x, highlighting why models with this final state 
tend to accommodate higher DM masses. 

Given the 0-step spectrum, determining the photon 
spectrum from an n-step cascade is particularly simple 
in the case of scalar mediators^] where the calculation 
essentially reduces to Lorentz-boosting the photon spec¬ 
trum up the ladder of particles appearing in the cascade. 
We review this calculation in Appendix |Bj As observed 
in J (>- , in the case of large mass hierarchies between the 
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steps in the cascade, the final photon spectrum can be 
simplified even further, as we now outline. 

Consider the zth step in the cascade, where the decay is 
4>i+i —> Let us define e* = 2mi/mi+i, and assume 

ei <C ljj Suppose the photon spectrum from decay of a 
single 4>i (and the subsequent cascade), in the rest frame 
of the 0, particle, is known and denoted by dN 1 /dxi_i. 
Then, in the presence of a large mass hierarchy, the decay 
of 4>i+i produces two highly relativistic 0j particles, each 
(in the rest frame of the 0i+i) carrying energy equal to 
TOi + i/2 = TOi/cj. The photon spectrum in the rest frame 
of the 02_j_i is then given by a Lorentz boost (see Ap¬ 
pendix fBI), and in the limit <C 1 takes the simple form 

m- 


dN ^ =2 f 1 ^_ 1 diV L+ 2 

dxi J x . Xi-i dxi-i 


(5) 


Here we have introduced the dimensionless variable Xi = 
2Ei/mi+i, where Ei is the photon energy in the 0i_|_i rest 
frame. Following this, once we know the 0-step spectrum 
we can iteratively derive the n-step result. The error 
introduced by this assumption is 0(ef), as we quantify 
in Appendix |B| 

Beyond simplifying calculations, the large hierarchy 
approximation is also convenient for the following two 
reasons. Firstly in this limit, we can specify the shape 
of the spectrum simply by the identity of the final state 
/, the value of ef, and finally the number of steps n. 
This is in contrast to the many possible parameters that 
could be present in a generic cascade. Secondly, as we 
will elaborate further in Sec. |IVj it is also possible to read 
off the results for a generic hierarchy once we know the 
small e.j result, making the assumption less limiting than 
it would initially appear. In particular in the limit when 
the masses become degenerate (e^ —► 1), the 00s are pro¬ 
duced at rest. When they subsequently decay, there is 
no boost to the 4>i+i rest frame, and so an n-step cascade 
effectively reduces to a hierarchical (n — l)-step cascade, 
except for the additional final state multiplicity. 

The Galactic frame is approximately the rest frame of 
the (cold) DM; consequently, to determine the measured 
photon spectrum, we need to calculate the photon spec¬ 
trum in the rest frame of the original DM particles. For 
an n-step cascade, this will involve n such convolutions, 
starting from the dN^/dxo 0-step spectrum, where the 
highest mass scale in the cascade will be mi =n = 2m x . 
Thus Xi =n = E n /m x , and the Galactic-frame photon 
spectrum will be dN 1 /dx n = m x dN 7 /dE n . Fig.^shows 
the resulting spectrum for a 0-6 step cascade in the case 
of final state taus with e T = 0.1. Each step in the cas¬ 
cade broadens out and softens the spectrum, and similar 
behaviour is seen for other final states. 


5 Note that the earlier-defined ej parameter does not function in 
exactly the same way as these e; parameters: ej fully parameter¬ 
izes the photon spectrum associated with production and decay 
of the SM particles, whereas the e, only describe Lorentz boosts. 
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We discuss the case of vector mediators in Sec. 
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In order to determine the favored parameter space, for 
a given choice of /, e/, and number of steps in the cascade 
n, we vary m x and an overall normalization parameter 
?7 (proportional to (<Tu)/m 2 , as we will see below) and 
compare the model to the data using the spectrum and 
covariance matrix of m■ In detail we calculate y 2 ac¬ 
cording to: 


X — ^ ) (A/),model A/) ; data) (A/j,model A/ydata) 5 

ij 

( 6 ) 

where 


model 

\ m x 

data = 

(E^ 
\ dE 


dN\ 
dx n ) 


i. model 


^ z,data 


(7) 

( 8 ) 


and both model and data are expressed in units of 
GeV/cm 2 /s/sr averaged over the region of interest. Here 
the C7i are elements of the inverse covariance matrix, 

1 J ^ 

which together with the data points are taken from m- 
By Eq.[2j the fitted normalization 77 is related to the DM 
mass and the J-factor by: 


(av) 


87T?7T 2 ?7 

J norm 


(9) 



10“ 3 10 -2 KT 1 10° 


X 


FIG. 2. An example photon spectrum from direct annihilation to 
taus (grey) and hierarchical cascades with n = ( 1 , 2 ,3,4,5, 6 ) steps, 
corresponding to (purple, blue, green, pink, orange, red) curves. 
The presence of each additional step in the cascade acts to broaden 
and soften the spectrum, and shift the peak to lower masses. All 
spectra are per annihilation. 


have enough mass to provide the rest masses of the decay 
products), setting a strict lower bound on the DM mass 
of: 


For consistency with the spectrum normalization of El 
the J-factor is averaged over the ROI \l\ < 20° and 2° < 
| 6 | < 20 °, so that: 

•/norm = [ dQJ ( l , b) / [ d<d 

J ROI J ROI (10) 

~ 2.06 18 x 10 23 GeV 2 cm” 5 . 


(Note that dfl = dldsinb, not dldcosb, since b measures 
the angle from the Galactic equator, not the north pole.) 

Self-Consistency Requirements: The procedure out¬ 
lined above treats m x as a free parameter that can be 
adjusted to modify the 0 -step spectrum; the fit only uses 
the shape of the spectrum provided by the 0 -step result 
and the boost of Eq. [5] However, there is an additional 
condition required for a cascade scenario to be physically 
self-consistent: the mass hierarchy between the DM mass 
and the particles produced in the final state must be suf¬ 
ficiently large to accommodate the specified number of 
steps. Equivalently, there is a hard upper limit on the 
number of steps allowed, for a given DM mass and final 
state. 

Recall that for an n-step cascade ending in a final state 
/, we defined e/ = 2 mf/mi , = 2 mi/m, 2 , e 2 = 270 , 2 / 1 x 13 

all the way up to e n = m n /m x . Combining these, the 
DM mass is given in terms of m/ and the e factors by: 


m x 


on rnj 


( 11 ) 


If the factors are allowed to float, we can still say that 
0 < < 1 in all cases (since each decaying particle must 


m x > 2 n m f /e f . 


( 12 ) 


In the remainder of this article we refer to this bound as 
a “self-consistency” condition or defining “kinematically 
allowed” masses. For consistency with the assumption of 
hierarchical decays (i.e. €i <C 1 ), the true bound on m x 
will in general be somewhat stronger than this conserva¬ 
tive estimate (although as we will discuss in Sec. 
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can become quite close to 1 before significantly modifying 
the fit relative to the —> 0 case). 


III. RESULTS WITH THE ASSUMPTION OF 
LARGE HIERARCHIES 

Here we present the results from the fits performed 
using the procedure outlined in the previous section. As¬ 
suming hierarchical cascades, we perform fits for four dif¬ 
ferent final states - electrons, muons, taus, and 6 -quarks 
- and fit over the photon energy range 0.5 GeV < E 1 < 
300 GeV0 Later in this section we discuss the effects 
of cutting out high energy data points, and how the fits 


By default, we omit the low energy data points with 0.3 GeV < 
< 0.5 GeV, as in this region the spectrum suffers larger 
uncertainties under variations of the background modeling, and 
the preferred value of the NFW 7 parameter is not robust m- 
We have confirmed that including these low-energy data points 
has little impact on our results. 
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FIG. 3. Contours of Ax 2 from the best-fit point (for a given 
step number n) corresponding to 1, 2 and 3 ct for final state 
fi’s, with fjj = 0.3. The purple, blue, green, pink, orange and 
red colors correspond to n = 1, 2, 3, 4, 5 and 6 steps in the 
cascades to final state //Is. Here we have fixed c ;( = 0.3 and fit 
over the range 0.5 GeV < E y < 300 GeV. 



FIG. 4. Contours of A\' 2 corresponding to 1, 2 and 3cr for 
n = 1 — 6 steps for e, //, t and b final states with £/ = 0.3. 
The fit is performed over the range 0.5 GeV < < 300 GeV. 

The best fit point of each step for all four final states follows 
a power law relation between m x and (/tv) . with index ~ 1.3. 
Only the darker regions are kinematically allowed. See text for 
details. 


would change if we only considered statistical uncertain¬ 
ties. 

In Fig. [3] we show a sample result, in which we plot 
Ay 2 1, 2 and 3cr contours in (m x , ( av )) space for 1-6 
step cascades ending in muons with = 0.3. The trend 
in the best fit point for each step is as expected. Recall 
the generic behavior illustrated in Fig. [2} each progressive 
step in the cascade acts to reduce the height of the peak 
and shift it to lower masses. Therefore higher steps in 
the cascades will be better fit by larger DM mass and 
cross-section as is indeed the case in Fig. [3] The larger 
cross-section results from an interplay of effects as can 
be seen from Eq. [9} an increased DM mass leads to a 
lower number density and hence a higher cross-section 
(scaling as to 2 ), but the increased power per annihilation 
implies a lower rj (adding a factor of m” 1 ), and finally the 
reduced height of the peak in the dimensionless spectrum 
for higher steps (as shown in Fig. [2j requires a larger 77 . 

In Fig. [ 4 ] we show the corresponding Ay 2 contours for 
electron, muon, tau, and 6-quark final states, again fixing 
Cf = 0.3. The best-fit mass and cross-section for each of 
the final states are empirically found to follow an approxi¬ 
mate power law with (av) oc m*' 3 . As discussed above we 
would expect (av) oc m x if the spectrum did not change 
in shape (simply being rescaled proportionally to m x to 
ensure energy conservation); the additional to 33 scaling 
factor reflects the change in shape of the spectrum. 

As discussed above, for a given DM mass and final- 
state fermion with mass rnf , there is an absolute upper 


limit on the number of steps allowed in a cascade, since 
every step corresponds to a change in mass scale of at 
least a factor of 2. In Fig. [IJ we show the contours if the 
limitation of Eq. [12] is ignored , since this conveys infor¬ 
mation on the mass scale and number of steps at which 
the broadness of the spectrum best matches the data; 
however, the mass values that violate this condition and 
so do not represent a self-consistent physical scenario are 
shown in lighter shading. This issue is relevant for the 
heavier final-state fermions, taus and 6-quarks, and par¬ 
ticularly acute for taus. Finally note that the irregu¬ 
lar shape of the contours for the one-step electrons and 
muons can be traced to the fact the 0-step FSR spectrum 
is both sharply peaked and has a kinematic edge, leading 
to a poor fit. 


In Fig. [H] we show the Ay 2 values between the best 
fit at a given step number n and the best fit overall, 
for each final state. We show results for both e/ =0.3 
and 0.1 in all cases, and also include e/ = 0.01 for elec¬ 
trons. As expected the results do not depend strongly 
on ej, especially in the case of taus, which is in accord 
with the results of Fig. [T] Note that the nominal overall 
best fit for the taus (n = 4) falls into the kinematically 
disallowed (inconsistent) region; n = 4 cannot be physi¬ 
cally accommodated within 3cr of its preferred DM mass. 
For this reason the results for taus and 6-quarks were re¬ 
run allowing only self-consistent scenarios (in the sense 
of Eq. 12); in these cases we obtain the results shown by 
the blue dotted curves in Fig. [5j We summarize the best 
fit results for e/ = 0.3 in Table |T] and the lcr range as 
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FIG. 5. Clockwise panels show the overall best fit for DM annihilating through an n-step cascade to electron, muon, 6 -quark and tau 
final states. The grey solid, dashed (and dotted) lines correspond to the Ax 2 between the best fit at that step, and the best fit for all n, 
for ef = 0.3, 0.1 (and 0.01) respectively. In the case of tau and 6 -quark final states, the blue dotted curves, denoted ‘physical,’ correspond 
to the case where only kinematically allowed (self-consistent) masses are considered as per the discussion in Sec. Ill] (we set ey = 0.3 for 
these curves). Note that in the case of taus, the “physical” best-fit points for 0 and 1 steps have the same x 2 as the best-fit points when 
“unphysical” scenarios are allowed, but as the overall best fit is different (with higher x 2 ) their Ax 2 is lower. The shaded bands correspond 
to the quality of fit. 0 -step results are not included for electrons and muons, as these fits are poor and have Ax 2 values well above the 
plotted y-axis. Electrons, muons and taus prefer longer 3-5 step cascades, whilst annihilations to 6 -quarks prefer shorter 0-2 step cascades. 
This is not surprising, since as has been already pointed out in the literature, 6 -quark final states are preferred for direct annihilations. 
Non-integer values of n can be associated with cascades containing steps with one or more large e^, as discussed in Sec. m 



FIG. 6. The blue, red, green and orange curves correspond to the 
overall best fit spectrum for e, //, r and 6 -quarks as determined 
from Fig. [5] Overlaid are the data points and systematic errors 
from nu. Note that due to correlations between energies, the best 
fit curves are not what would be naively expected if only statistical 
errors were present. 

determined from Fig. [5] on these parameters in Table [TT| 
In Fig. [g] we show the overall best fit spectrum for elec¬ 
tron, muons, taus, and 6 -quarks with e/ = 0.3. Although 
the spectra for direct annihilation to these final states 
are quite different, after introducing the freedom to have 


Final State 

n-step 

m x (GeV) 

<tv (cm 3 /sec) 

x 2 

e 

5 

67.2 

2.9 x 10 -24 

26.82 


4 

53.0 

9.9 x 10~ 25 

26.94 

^unphysical 

4 

59.4 

4.6 x 10 -26 

24.13 

^physical 

2 

24.1 

1.4 x 10" 26 

25.59 

b 

2 

91.2 

3.9 x 10 -26 

22.42 


TABLE I. Best fit to DM annihilations to various final states 
with tf = 0.3. For the case of taus we show a best fit point if 
we include kinematically disallowed masses (unphysical) and 
also if we restrict ourselves to physical masses as discussed in 
Sec. |TT] Fits were performed over 20 degrees of freedom. 


multi-step cascades, a similar best fit spectrum is picked 
out in each case. To expand on this, we can compare 
the various 0 -step spectra - as displayed in Fig. [T] - to 
the result of a hierarchical n-step cascade that ends in 
cf >i —> 77 . This comparison is shown in Fig. [7] The spec¬ 
trum of photons from this process is just a 6 -function in 
the 4>\ rest frame, and is in a sense the simplest possi¬ 
ble photon spectrum. We find that the photon spectrum 
from direct annihilation to electrons is similar to that ob¬ 
tained by a 2-3 step cascade terminating in cj>i —► 77 ; for 
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Final State 

n-step 

m x (GeV) 

ov (cm 3 /sec) 

e 

3-6 

28-107 

10“ 24 °-10“ 23 - 3 

M 

2-5 

22-89 

io- 245 -io- 23 - 7 

Funphysical 

3-5 

37-94 

io- 256 -io- 251 

^physical 

2 

24.1 

IQ- 25 - 8 

6 

0-3 

40-150 

j^q — 25.8 j^Q—25.2 


TABLE II. Range of parameters within lo of the best fit step 
for e/ = 0.3 for electrons, muons, taus and 6-quarks. As in 
Table [I] we show both physical and unphysical tau results. 

muons and taus the closest match is a 3-4 step cascade; 
and for 6 -quarks 6-7. Of course these correspondences are 
not exact - for example, the 6 -quark spectrum is more 
complex than just applying Eq. [5] to a 6 -function - but 
they allow us to regard these 0 -step spectra as arising ap¬ 
proximately from a common ( 6 -function) spectrum con¬ 
volved with differing numbers of cascade steps. We can 
then intuit how many additional steps are required in 
each case, to bring the spectra to a similar shape. Com¬ 
bining these numbers with the preferred number of steps 
seen in Table [I] we find the GCE prefers a spectrum that 
can be roughly modeled as a 6 -function occurring at the 
endpoint of 7-9 cascade decays. In this sense it seems 
fits to the GCE prefer a cascade with a large number of 
steps, and that these can occur in the SM or dark sector. 

Likewise, this general picture can approximately de¬ 
scribe showers in the dark sector [30]. Such showers will 
effectively contain decay cascades of different lengths, but 
we find that the spectrum of [3D] can be well described by 
a 6 -function </>i —► 77 broadened by ~ 3 decay steps. The 
best-fit scenario found in that paper corresponds to a DM 
mass of ~ 10 GeV; this is consistent with the preferred 
mass for our 1 -step electron case, which also corresponds 
to a 6 -function at the endpoint of a ~ 3-step cascade. 
A better fit to the data might therefore be obtained by 
combining such dark showering with a short dark-sector 
cascade. In Sec. |IV| we will return to this point, and 
discuss the sense in which our results may be used to 
estimate the parameter space for dark shower models. 


A. Different Final States 

A few comments about the various final states are in 
order. 

Electrons: The photon spectrum from direct annihi¬ 
lations XX ~> e + e _ is sharply peaked. This tends to 
produce a worse fit to the GCE. As such we need several 
steps in the cascade in order to broaden the spectrum 
sufficiently to allow for a parameter space where a sig¬ 
nificantly improved fit is possible, and this is shown by 
the substantial decrease in the quality of fit at low n in 
Fig .0 It should be noted that any model for the GCE 
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FIG. 7. The 0-step spectra for e, /i, r and 6-quarks with £/ = 0.3 
are shown as the blue, red, green and orange curves. The dashed 
curves show the spectrum of a hierarchical n-step cascade that ends 
in 0i —» 77 (a 5-function in the 0i rest frame) for n = 1 — 7, with 
lighter curves corresponding to progressively longer cascades. In 
order to compare the shape of the spectra we have magnified the 
0-step spectra by a factor of 470, 190, 6.2 and 3.1 for e, p, r and 
6 -quarks respectively. We see the electron spectrum is closest to a 
2-3 step cascade ending in a 5-function, muons and taus are closest 
to a 3-4 step cascade, whilst 6-quarks most resemble 6-7. 


with direct annihilation into electrons will likely be in se¬ 
vere tension with the data from AMS [3SJ. This tension 
is likely to persist for at least the n = 1 cascade, and pos¬ 
sibly higher steps as well [HDj. As we go to higher-step 
cascades the spectrum broadens and the AMS bounds 
are expected to weaken, but the exact bounds should be 
worked out for any cascade scenario with a branching 
fraction to electrons. For the purposes of this work, we 
use the electron case as an example of a sharply peaked 
photon spectrum to demonstrate the impact of the spec¬ 
tral broadening, not necessarily as a realistic explanation 
for the excess. Similarly, constraints on DM annihilation 
from the cosmic microwave background (CMB) (DD] are 
likely to rule out both the electron and muon favored re¬ 
gions shown in Fig. [Tj while leaving the 6 and tau regions 
largely unconstrained. The figure of merit for CMB con¬ 
straints is ( av)/m x [TOFID] , up to an Oil) factor which is 
channel- and spectrum-dependent [13101]. As discussed 
above, for the best-fit regions (for hierarchical decays), 
this quantity scales as ~ m^ 3 as the number of steps in¬ 
creases; thus, we expect the constraint to become slightly 
stronger for longer cascades. 

Muons: In Fig. [5] we see that the muon final state spec¬ 
trum has the same qualitative behavior as the electrons, 
and will be subject to similar constraints. This is un¬ 
surprising as the muon spectrum is quite similar to that 
from electrons, albeit with a less pronounced peak (see 
Fig. [TJ. 

Taus: As with other leptonic final states, taus also 
prefer multi-step cascades for the best fit. Note that the 
best fit point at 4 steps is in fact kinematically disallowed 
(inconsistent) as can be seen in Fig. [I] and as discussed 
















FIG. 8 . The 3cr contours for 1-6 step cascade annihilations to 
final state electrons with e e = 0.1. Red contours correspond to 
fitting over the entire energy range 0.5 GeV < F 7 < 300 GeV with 
the full covariance matrix of ED- Orange contours correspond to 
fitting with a cut on high energies 2£ 7 < 10 GeV. Green contours 
correspond to a fit over the full energy range but with only the 
statistical errors of [Ill- 

in Sec. [XT] However, the best fit point after imposing the 
consistency condition, at 2 steps, is still a better fit than 
the high-step cases with electron and muon final states. 

b-quarks: DM annihilation to 6-quarks is the preferred 
channel for direct annihilation identified in mm , where 
it already provides a good fit. Accordingly there is no 
need to broaden the spectrum with a large number of 
cascades - however, as we will discuss in Sec. [Vj even a 
short cascade can greatly alleviate constraints from col¬ 
liders and direct searches (see also f27] i28j and references 
therein). A cascade with several steps can still give an 
equally good or slightly better fit, and of course accom¬ 
modates higher masses than for the case of direct anni¬ 
hilation. However, since the spectrum is already fairly 
broad, adding too many additional steps makes the fit 
worse, as shown in Fig. [5] Accordingly, the DM mass can¬ 
not be pushed far above 100 GeV without significantly 
worsening the fit, at least in the context of hierarchical 
cascades. 


B. Sensitivity to Systematics and Energy Cuts 

In the results presented above we have fit the data 
of m to the photon spectrum from DM annihilations 
through multi-step cascades to various final states. The 
fit was performed over the energy range 0.5 GeV < E 1 < 
300 GeV. There is some evidence that the emission de¬ 
tected above 10 GeV may not share the same spatial 


profile as the main excess, suggesting a possible indepen¬ 
dent origin (for example, these high-energy data appear 
to prefer a morphology centered at negative l and with 
a shallow spatial slope HU), so we also test the impact 
of omitting the data above 10 GeV. Finally, we explore 
the impact of including only the statistical uncertainties 
of m, omitting systematic errors, to test the degree to 
which the constraints could improve with reduction in 
the systematic uncertainties. 

We display the results of this study in Fig. [8][9j for the 
case of n-step cascade annihilations to final state elec¬ 
trons with e e = 0.1. Annihilations to other final states 
generically display the same behavior as the energy range 
and error estimates are varied. Cutting out the high en¬ 
ergy data points generically shifts the fit to prefer lower 
masses and narrower spectra, and therefore corresponds 
to cascades with fewer steps - resembling a 6-function at 
the endpoint of a 5-7 step cascade, rather than a 7-9 step 
cascade. At a fixed number of steps, the main impact 
of omitting the high-energy data points is to raise the 
preferred cross-section and shrink the contours. Under¬ 
standing the high-energy data will thus be important in 
distinguishing quantitative models for the GeV excess. 

Fitting over statistical errors increases the actual \ 2 
values, and the rate at which \ 2 increases away from its 
minimum (as expected), as demonstrated by the shrink¬ 
ing green contours of Fig. [8j The overall preferred step 
in the cascade however is not dramatically affected, only 
changing by 0-1 steps, as shown in the top panels of Fig. [9] 
- we display the corresponding best fit spectra in the bot¬ 
tom panels. At a fixed number of steps, the preferred 
cross-section increases, becoming more similar to what 
we find when omitting the high energy points. 


IV. INTERPRETATION FOR GENERAL 
CASCADES 

A. Relaxing the Assumption of Large Hierarchies 

The results displayed in the previous section were ob¬ 
tained assuming large mass hierarchies between each cas¬ 
cade step. It is possible to recast these results to gain in¬ 
sight into the case of general £j values. To see this, con¬ 
sider the decay <j>i+i — > 4>i4>i- As previously discussed, 
in the limit when two mass scales become degenerate 
(ej -A 1), an ?r-step cascade effectively reduces to an 
(n — l)-step cascade, except for the additional final state 
multiplicity. Thus adding a degenerate step to a cascade 
is much simpler than adding one with a large hierarchy: 
we need only multiply the spectrum by two to account 
for the increased multiplicity, and halve the photon en¬ 
ergy scale to account for the initial energy being spread 
between twice as many particles. (For completeness, we 
check analytically that the limit of e* —> 1 has this be¬ 
havior in Appendix [b|) 

In light of this, an n-step cascade with one degenerate 
step and an (n— l)-step hierarchical cascade must provide 
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FIG. 9. Top Panels: We show the impact on the preferred number of steps when changing the energy range and error types considered. 
Each curve is for final state electrons with e e = 0.1. The left figure shows the use of systematic errors over the full and a restricted energy 
range (E 7 < 10 GeV) in red and orange respectively. The right figure is the equivalent for statistical errors, with the full energy range 
shown in green and the restricted in blue. Bottom panels: Here the best fit curves as determined from the top panels are shown with 
the appropriate data and errors from m overlaid, for the example case of the electron final state. The left panel shows the results for 
systematic errors, where the best fit point was n = 5 for the full range (red curve) and n = 3 for the restricted range (orange curve). The 
right panel shows the equivalent for statistical errors, where for the full range the n = 6 curve is shown in green and for the restricted 
range the n = 2 curve is in blue. 


equally good fits to the GCE, with the former preferring 
twice the annihilation cross-section and DM mass rela¬ 
tive to the latter. The increased DM mass results from 
the halving of the energy scale, whilst to understand the 
cross-section we look back to Eq. [9] adding the degener¬ 
ate step doubles the photon multiplicity, which halves tj 
to compensate, but the doubling of the DM mass means 
overall the cross-section is increased by a factor of two. 
As such the results in Fig. [4] can be readily extended for 
additional degenerate steps. For each additional degen¬ 
erate step on top of an initial hierarchical cascade (the 
degenerate step may occur anywhere in the cascade), the 
shape of the x 2 contours remains the same, but shifted 
upward by a factor of two in mass and cross-section. 
With a sufficiently large number of degenerate decays, 
the DM mass required to fit the GCE could be made ar¬ 
bitrarily high, although this would seem to require con¬ 
siderable fine-tuning. (A natural scenario in which one 
degenerate step arises due to a symmetry is discussed in 

0510 

Cascades with general values of e, in turn interpolate 
between the two simpler cases already considered, with 
small and large e,. We give the general convolution for¬ 


mula in Appendix [Bj and an example of how spectra 
evolve as a single shifts from 0 to 1 is shown in Fig. [lO] 
This interpolation provides an alternate interpretation 
for Fig. [5] the n on the z-axis of these plots can be 
thought of as representing the number of steps with a 
large hierarchy, rather than the total number of steps. 
If one of these steps becomes degenerate (while holding 
the total number of steps fixed), as previously discussed, 
we will move from n to n — 1 steps in terms of the spec¬ 
tral shape and hence quality of fit. Intermediate Ci val¬ 
ues will interpolate smoothly between these two cases. 
Thus for any arbitrary collection of hierarchical and de¬ 
generate steps, the quality of the fit and the location of 
the best-fit region in m x — (av) parameter space can al¬ 
ready be estimated from Figs. [l][5j A concrete example 
of the transition in preferred DM mass and cross-section 
is shown in Fig. El which corresponds to the variation 
of the spectrum shown in Fig. [lOj The curve plotted out 
by the best fit point for intermediate values of e is not 
a straight line between the two extreme values, but does 
not deviate far from this. Similar behavior was seen for 
other final states and choice of degenerate step. 

At a fixed DM mass, the perturbation to the = 0 
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FIG. 10. The transition of the spectra between € 2=0 and €2 = 1, 
calculated using Eq. |B5| The example case is a 2-step cascade with 
final state taus and e r = 0.1. The dark blue is for e = 0 and is what 
would result from the large hierarchies approximation. The e = 1 
case shown in light blue corresponds to a completely degenerate 
spectrum, and as such is equivalent to a shifted 1-step curve. In 
between these two, we show intermediate e values as the dashed 
curves, specifically e = {0.3,0.5,0.7,0.9,0.99}. Note the rate of 
transition between the two cases is in keeping with the error in the 
large hierarchies case being of order 0(e|). 



FIG. 11. The transition of the best fit point and lcr contours 
between €2=0 and €2 = 1, calculated using Eq. |B5| The example 
case is a 2-step cascade with final state taus and e T = 0.1. The 
transition is between the e = 0 case in dark blue and e = 1 in light 
blue. The dashed curves map out the transition with intermediate 
values, specifically e = {0.3,0.5,0.7,0.9,0.99}. 


photon spectrum evolves roughly as ef as varies from 
0 to 1 (as discussed in Appendix [B]); this behavior is 
shown in Fig. [l0| where the £2 = 0.3 spectrum is almost 
indistinguishable from the £2 = 0 spectrum, and £2 = 0.5, 
£2 = 0.7 and £ 2 = 0.9 give spectra intermediate between 
the £2 = 0 and £2 = 1 cases. The perturbation to the 
best-fit y 2 will tend to increase even more slowly than 
ef, in the case where e^ = 0 is a better fit than e, = 1 , 
since the DM mass and cross-section can float to absorb 
changes in the spectrum and reduce the increase in y 2 . 
In all examples tested the best-fit y 2 remains essentially 
unchanged from the e, = 0 case out to e t = 0.7. 

In general a cascade with n total steps, rid of which are 
degenerate (rid values of e^ —> 1 ) will have the same spec¬ 
trum as a cascade with (n — rid) hierarchical steps with 
a factor of 2 nd enhancement in mass and cross-section. 
This is illustrated in Fig. [T2] for the case of decays to fi¬ 
nal state t’s with 1-6 total cascade steps. Relaxing the 
assumption of large hierarchies therefore results in a pre¬ 
ferred triangular slice of parameter space, bounded by 
curves with (av) oc m x and (av) oc m 3 ' 3 . We can now 
understand the results of Fig. [5] as mapping out the vari¬ 
ation in y 2 when moving between classes of scenarios, 
each defined by a fixed number of hierarchical steps but 
containing scenarios with varying numbers of degenerate 
steps (each of these classes is represented by a line in 
Fig. [i~2|. Note also that the kinematic constraint Eq. 11 
acts on classes rather than individual scenarios (since 
adding a degenerate step doubles the DM mass but in¬ 
creases the number of steps by 1 , strengthening the con¬ 
straint on DM mass by a factor of 2); if one scenario is 
disallowed the entire class is disallowed. 


Fig. 13 summarizes our combined results. There, 


the top panels display the regions mapped out in the 
(av) — m x plane by the best fit points involving 1-6 steps 
(either hierarchical or degenerate) cascades to final state 
electrons, muons, taus and 6 -quarks. In the bottom pan¬ 
els, we indicate which hierarchical step and final state 
yield the best fit, and the comparative quality of fit for 
other combinations. We show all these results for fits 
over the full (left panels) and restricted (right panels) 
energy ranges. Additionally as shown in the top panels, 
electrons (taus) and muons ( 6 -quarks) have some degree 
of overlap, especially once degenerate steps are included. 
The overlap of these regions is reduced when the high en¬ 
ergy data points are excluded, as is clear by comparing 
the right and left panels. 

The positions of the triangular regions in Fig. [T3] 
largely reflect the differing branching ratios to photons 
(rather than other stable SM particles) for the different 
final states. For each of the direct annihilation (0-step) 
spectra, we can compute a factor fc, defined as the to¬ 
tal energy in photons (per annihilation) as a fraction of 
m, = 2m x . For example, direct annihilation/decay to 
77 would have k = 1. For the final states we consider, 
we find k = 3.0 x 10 -3 , 7.0 x 10 —3 , 0.14 and 0.26 for 
electrons, muons, taus and 6 -quarks respectively. Final 
states with smaller k will naturally require higher cross- 
sections in order to fit the signal. In Fig. [If] we show the 
results of Fig. 13 replotted in terms of k(av) and m x : 
we see that once this factor is taken into account, all 
channels pick out essentially the same triangular region 
of parameter space, bounded by curves with k(av) oc m x 
and k(av) oc m 3 ' 3 . 

Incorporating dark showers: This concordance be- 
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FIG. 12. The purple, blue, green, pink, orange and red points cor¬ 
respond to the best fit ( m x , crv) point for a total number of cascade 
steps (degenerate + hierarchical) n = 1, 2, 3, 4, 5, 6 respectively; 
for annihilations to final state taus with e T = 0.3. Points living 
on the same line have the same number of hierarchical steps and 
therefore result in equally good fits to the data. Points of the same 
color, but with progressively greater values of ( m x , crv ), correspond 
to successively replacing hierarchical steps with degenerate steps, 
holding the number of total steps fixed. For the above case of taus 
only the one and two step hierarchical cascades are kinematically 
allowed as indicated in Fig. |4] (note that the kinematic constraint 
applies to lines as a whole, not individual points; see text), thus 
only points living on the solid lines are allowed as these lines corre¬ 
spond to cascades with one and two hierarchical steps respectively. 


tween the different final states suggests that dark shower 
models may be expected to also inhabit this region. For 
instance, the authors of j30] find a preferred cross-section 
of 8x 10~ 27 cm 3 /s for their SU(2)y model, with a roughly 
35% branching ratio into stable dark sector baryons (with 
other decay channels ending in photons), and a preferred 
mass of ~ 10 GeV. At first glance this suggests a some¬ 
what higher value for k(av ) than the lower tip of the 
triangular region identified in Fig. 14 However, [30] fits 
to a different spectrum for the GCE excess (taken from 
[10]). without a systematic uncertainty estimate, and as¬ 
sumes a lower local DM density (0.3 GeV/cm 3 rather 
than 0.4 GeV/cm 3 ) Q In our analysis, omitting system¬ 
atic errors (or removing high-energy data points) raises 
the preferred cross-section by a factor of ~ 2 (Fig. [8]), 
and likewise lowering the local DM density from 0.4 to 
0.3 GeV/cm 3 would raise the required cross-section by a 
factor of ~ 2; the lower tip of our triangular region would 
then reside at m x ~ 10 GeV and k{av) ~ 3 x 10 -27 
cm 3 /s, which seems roughly consistent with [3Dj . 


7 Private communication, Dean Robinson. 


B. Models with Vector Mediators 

Thus far we have considered models of multi-step 
cascades through scalar mediators. However models 
in which the hidden sector mediators include vector, 
fermion or pseudo-scalar particles are at least as equally 
well motivated (e.g. |26] or the dark shower example dis¬ 
cussed above [3D]). In the case of vector or fermionic me¬ 
diators the simple recursion formula Eq. [5] will in general 
no longer hold, since the photon spectrum from the decay 
of mediators with spin need not in general be isotropic. 
The standard recursion formula will also break down if a 
decay is more than two-body, or if the decay is two-body 
but the decay products have different masses (although 
if the decay is strongly hierarchical the impact will be 
tiny), since these possibilities modify the Lorentz boost 
from the fa frame to the fa+i frame. Note this is differ¬ 
ent to having several possible decay chains with different 
branching ratios; in this case our analysis does apply, and 
the final spectrum will simply be a linear combination of 
the spectra produced by the different decay chains. 

Anisotropy of the photon spectrum is not in itself a 
sufficient condition for the recursion formula to break 
down. To modify the recursion, for some step i, the dif¬ 
ferential decay rate of fa must be a function of the angle 
9 between (1) the momenta of the decay products in the 
c j>i rest frame and (2) the boost direction from the fa 
rest frame to the <pi-\ rest frame. (Here we use fa to 
denote arbitrary mediators, independent of their spin.) 
Since the decays in the fa rest frame do not “know” 
about the fa+± frame, this sort of correlation is only pos¬ 
sible if (1) the direction of the spin/polarization vector of 
the fa in its rest frame depends on the momentum with 
which it was produced in the rest frame, and (2) 
the spectrum of the decay products of fa is a function of 
the angle between their momentum and the rest-frame 
spin/polarization vector of fa. If only one of the two 
applies, averaging over the spin/polarization of fa will 
leave no 0-dependence. However, both these properties 
will generically hold if fa is a vector: typically the decay 
of fa- 1 will prefer either longitudinally or transversely 
polarized vectors 4>i , which will in turn decay with differ¬ 
ent angular distributions. 

Let us consider the potential impact of such a 9- 
dependence. For illustrative purposes, let us suppose 
that the photons produced in the decays of fa (whether 
directly or by subsequent decays of the fermions) have es¬ 
sentially the same energy spectrum as in the pure-scalar 
case, in the rest frame of the fa. This assumption might 
fail if the spin of fa affects the correlations (if any) be¬ 
tween the fermion spins, fermion momenta and photon 
momenta, but by making it we can isolate the impact of 
angular dependence in a single step of the cascade. 

Consider a one step cascade XX ► fa fa, fa —► //, 
where fa is a vector boson. Suppose the full spec¬ 
trum of photons in the fa rest frame can be written as 
^ = /o (yo) dN/dx 0 , where y 0 = cos 9 0 and dN/dx 0 is 
the spectrum for the scalar mediator case fa = 1. Then 
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FIG. 13. Combined results of fits with ey = 0.3 over the full energy range (left) or with a restriction < 10 GeV (right). Top panels: 
Best fit ( m x ,(7v ) for a cascade with 1-6 total (degenerate + hierarchical) steps ending in electrons, muons, taus of 6-quarks. Points on 
the sa me l ine have the same number of hierarchical steps and therefore result in equally good fits to the data, following the discussion 
in Sec. |IV| Points of the same color, but with sequentially greater values of ( m x ,(7V ), correspond to progressively replacing hierarchical 
steps with degenerate steps, holding the total number of steps fixed. The color of the lines indicate goodness of fit and only solid lines 
are kinematically allowed (as explained in see Sec. |h|). Bottom panels: Show the overall best fit for DM annihilation through an n-step 
hierarchical cascade to electron, muon, tau and 6-quark final states. The curves show the Ax 2 of the best fit at that step and final state, 
as compared with best fit over all steps and final states. No restriction to physical kinematics is imposed, but where restrictions would 
apply can be inferred from the top panels. The shaded bands correspond to the quality of fit. For fits over the full energy range a fairly 
short cascade terminating in decay to 6-quarks gives the preferred spectrum, whilst over the restricted energy range each final state can 
potentially provide approximately equally good fits. 


the now familiar formula for the energy spectrum in the 
XX center of mass frame is: 



where we calculated the yo integral assuming ei <C 1. 
Again we could extend this expression to an n-step cas¬ 
cade using the same formalism as in Appendix [B] The 
angular dependence at each step will in general be dif¬ 
ferent depending on the model; we can parameterize this 
by specifying different functions fi (yi) at each step. In 


the limit of small we find: 


dN 1 _ 2 f 1 dXi-1 
dxi J Xi Xi—\ 

(14) 

A detailed study of the impact of vector or fermionic 
mediators is beyond the scope of this paper; we leave it 
to future work. However, we can work out an explicit 
example motivated by the case where at the end of the 
cascade, a scalar/pseudoscalar resonance decays to two 
vectors which subsequently each decay into two fermions. 
This scenario has been studied in the context of Higgs 
decays J4B], furnishing results for a general resonance X 
decaying to two identical vectors VV, which each in turn 
subsequently decay to //. (In our notation, the V here 
would correspond to </>i and X to </> 2 .) The differential 
decay rate to fermions in this case is a linear combination 
of terms proportional to sin 2 9, 1+cos 2 9 and cos 9 (where 


2 Xi 

Xi -1 


- 1 


dN 

dxi-i 


+ <3(e 2 ). 
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FIG. 14. Colored points indicate the best fits for different numbers 
of hierarchical and degenerate cascade steps, and different final 
states, as in Fig. |l3| However, here we rescale the cross-section 
by the fraction of power into photons k for each final state (3.0 X 
10 -3 , 7.0 x 10 -3 , 0.14 and 0.26 for electrons, muons, taus and b- 
quarks respectively). All final states then pick out the same region 
of ( m x ,kcrv ) parameter space. The dashed lines indicate curves 
with k(av) oc m x and k(crv) oc ra*- 3 , chosen to originate from the 
lowest-mass point studied; these curves approximately bound the 
full parameter space of interest (see text). 


9 is the angle defined above and in Appendix [b]), with co¬ 
efficients depending on the axial and vector couplings of 
the fermions to the V, and the parity of the initial state 
X Hi]. In hierarchical decays of a scalar or pseudoscalar 
resonance to VV, where V has vector (rather than axial 
vector) couplings to //, the dominant angular depen¬ 
dence is either 1 + cos 2 9 or sin 2 9. For these specific (but 
common) angular dependences in the <j>i decay, we show 


the resulting changes to the photon spectrum in Fig. 15 


The impact is modest, and so we expect our qualitative 
results should hold for more general cascades. 


V. SIGNALS AND CONSTRAINTS 

While we have remained agnostic regarding the choice 
of an actual model, we point out that any model with 
new light states in a dark sector that explains the GCE 
must also be consistent with the following experimental 
constraints: 

• Direct Detection: The coupling controlling ctdd 
must not be so large as to be in conflict with bounds 
from DM direct detection experiments m- 

• Big Bang Nucleosynthesis (BBN): New light states 
must decay fast enough such that they do not spoil 
the predictions of BBN. 
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x 

FIG. 15. Spectrum for a 1-3 step cascade with a vector mediator 
in the final step of the cascade <f >2 —> VV, V —> //. We consider 
three separate cases: f(9) = 1, (3/4)(l 4- cos 2 0), and (3/2) sin 2 9. 
The first of these is equivalent to a cascade with only intermedi¬ 
ate scalars (and hence isotropic decays), the others correspond to 
common angular dependences (see text). 


• Collider constraints. 

These experimental constraints on a multi-step cascade 
will be very similar to those on a one-step cascade, with 
the key parameter being the coupling of the dark sector 
to the SM in both cases. 

The simplest models that explain the GCE by direct 
DM annihilations to SM states are generally in conflict 
with direct detection bounds: the same coupling that 
must be small enough to avoid the LUX bound m , 
must also be large enough to explain the GCE with a 
thermal WIMP (note however that this conclusion is not 
inevitable; there are effective DM-SM couplings and sim¬ 
plified models that generically evade the bounds, e.g. 
[2T1 [22] ). As pointed out in [261428] , the addition of a dark 
sector with a single mediator allows for an explanation 
of the GCE while alleviating direct detection constraints. 
The reason is straightforward: any direct detection signal 
will be controlled by the coupling of the mediator to the 
SM, whereas the annihilation rate is independent of this 
quantity, so the two can be tuned largely independently. 
We make this point more explicit in Appendix [C] Ex¬ 
actly the same property holds in models with expanded 
cascades, where the direct detection signal is controlled 
by the coupling between the dark sector and the SM; 
indeed, the direct detection signal may be suppressed 
even further if the coupling between the DM and the 
SM requires multiple mediators. If the couplings within 
the dark sector are not highly suppressed, decays within 
the dark sector should in general proceed promptly (on 
timescales C 1 s), and so the constraint from BBN will 
primarily limit the coupling of the final mediator in the 
cascade to the SM. Accordingly, since it has been shown 
that for one-step cascades the constraint from BBN can 
be consistent with a null signal in direct detection ex¬ 
periments m, the same should hold true for multi-step 
cascades (since in the multi-step case, the final step con- 
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trols the coupling to the SM and hence provides the only 
relevant parameter for both BBN and direct detection). 
Collider bounds and limits from invisible decays of SM 
particles are also controlled by this final coupling, so can 
accordingly be dialled down in the same way as for one- 
step cascades, consistent with BBN bounds on the final 
coupling m- A complex dark sector with multiple me¬ 
diators could potentially give rise to interesting collider 
signatures (e.g. .'17 IB HU j • but a detailed discussion is 
beyond the scope of this work. 


VI. CONCLUSION 


spectrum associated with the excess; we do not take a po¬ 
sition on this question, but note that its resolution may 
affect the range of dark-sector models that can provide 
viable explanations of the excess. 

In this work we assumed that the directions of decay 
products in the rest frame of their progenitor are uncor¬ 
related with the direction of the Lorentz boost to the 
rest frame of the previous progenitor particle in the se¬ 
quence. Whilst always true for scalars, this may not hold 
for vector and fermionic mediators. We leave a more de¬ 
tailed discussion of concrete multi-step cascade models 
exploring these issues for future work. 


We have laid out a general framework for characteriz¬ 
ing the photon spectrum from multi-step decays within 
a secluded dark sector terminating in a decay to SM par¬ 
ticles, and explored the ability of such a framework to 
produce the GeV gamma-ray excess observed in the cen¬ 
tral Milky Way. 

For any given SM final state, allowing multi-step de¬ 
cays expands the preferred region of m x — (av) to a trian¬ 
gular region of parameter space, probed by cascades with 
different numbers of degenerate and hierarchical decays 
(where the decay products are slow-moving or relativis¬ 
tic, respectively), and bounded by curves with (av) oc m x 
and (av) oc to*’ 3 . Decays to different Standard model fi¬ 
nal steps correspond to different triangular regions in pa¬ 
rameters space as shown in Fig. |13[ Large numbers of de¬ 
generate decays can raise the mass scale for the DM with¬ 
out bound, albeit at the cost of requiring a cross-section 
much higher than the thermal relic value and some de¬ 
gree of fine-tuning. Hierarchical decays broaden the pho¬ 
ton spectrum, permitting a better fit to the data for SM 
final states that produce a sharply peaked photon spec¬ 
trum; however, more than 4-5 hierarchical decays begin 
to reduce the quality of the fit even if the initial spectrum 
is very sharply peaked. In the absence of degenerate de¬ 
cays, the preferred mass range for the DM can then be 
constrained, and is consistently ~ 20—150 GeV across all 
channels; the corresponding cross-sections are close to the 
thermal relic value for tau and 6 -quark final states, and 
1-2 orders of magnitude higher for e and /i final states. 
Regardless of the final state, with the additional freedom 
of hierarchical decays the preferred spectrum tends to a 
similar shape, which can be approximated as the result 
of a cascade of 7-9 hierarchical decays terminating in a 
two-body 77 decay. We find that the best overall fits are 
still attained by DM annihilating to 6 -quarks (or other 
hadronic channels) with 0-2 hierarchical steps. 

Our preferred (av) — m x regions are fairly insensitive 
to the details of the uncertainty analysis or the range 
of data points included. However, omitting high-energy 
data (above 10 GeV) substantially reduces the preferred 
number of hierarchical decay steps (from 4-5 to 2) for 
channels where the photon spectrum from direct annihi¬ 
lation is sharply peaked. There is currently disagreement 
between different analyses as to the high-energy photon 
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Appendix A: O-step Spectra 


In order to calculate the photon spectrum, it is more 
straightforward to first determine the density of states 
according to: 


annihilations : 


decays : 


1 diV 7 
iV 7 dE-y 

1 dNy 
iV 7 dE-y 


1 d(av) 
(av) dE-y 

i dr 

r dEry 


(Al) 


from which the spectrum can be easily backed out. Note 
that as pointed out in [50], if the cascade begins with 
a decay x 4>n(f>ni we will obtain an identical photon 
spectrum to the annihilation scenario, except the initial 
DM particle will be twice as heavy. This is the sense in 
which our results are readily transferred to the case of 
decaying DM. The key difference for the decaying case 
is the spatial morphology of the signal will generically 
require a line of sight integral over the DM density, rather 
than density squared as appears in the J-factor in Eq. [3] 
The observed spatial morphology of the GCE appears 
to disfavour decaying scenarios, which is why we do not 
mention them further here, although see m for a novel 
decay scenario that is distributed like density squared. 

The result of Eq. m is that in some circumstances it 
is possible to calculate various step cascades analytically. 
This approach is shown for several cases in 150]. Yet 
in many cases - most notably those involving hadronic 
processes in their final states - analytic calculations are 
not feasible. For the present work we used a combination 
of analytic and numeric results depending on the final 
state employed. The details for each case is outlined 
below. 
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1. Annihilations to e + e 

The only contribution to the photon spectrum arises 
from FSR via the decay (f>i —> e + e~ r y. The spectrum in 
this case can be calculated analytically using Eq. |A1[ 
which was done in [36j for the generic case of (f>\ —> 
f + f~ 7 . As pointed out there, when using the simple 
convolution formula Eq. [5j consistency requires throw¬ 
ing away terms 0(e 2 ) and higher, where e/ = 2mf/rrii. 
Doing so they obtained the following expression for the 
spectrum that we include for completeness: 

dN RSR _ «em 1 + (1 — x 0 ) 2 

dx 0 7T Xq 

(A2) 

Note the In term will dominate for small €f, and the — 1 is 
simply included to ensure consistency with the large hier¬ 
archies approximation. We confirmed that this spectrum 
is in agreement with the output from Pythia8 in the case 
of final state electrons. From here, by repeated use of the 
convolution formula it is possible to obtain completely 
analytic formula for the ro-step cascade, which were used 
in our fits. For example, the first two steps are shown in 

ESj. 

2. Annihilations to fi + 


hadronize (dominantly to pions) which will result in large 
contributions to the photon spectrum. We simulated this 
final state in Pythia8 to generate an initial spectrum, to 
which we could then apply the convolution formula. 

4. Annihilations to bb 

Much like for taus, in the case of final state 6 -quarks 
FSR is a subdominant contribution, and instead the spec¬ 
trum is largely determined by hadronic processes. As 
such we again utilize Pythia8 to obtain the initial spec¬ 
trum. 


Appendix B: Kinematics of a Multi-step Cascade 

As already emphasized the utility of the small = 
- or large hierarchies - approximation is three¬ 
fold: 

1. It simplifies calculations in that we can use Eq. [5] 
rather than the general formula we display below; 

2. More importantly it allows us to describe a cascade 
using just the identity of the final state /, the value 
of e/, and the number of steps n, in contrast to the 
many possible parameters of the generic case; 




For final state muons, in addition to FSR, as pointed 
out in |36j the radiative decay of the muon yL —> eVeV^ 
will meaningfully contribute to the photon spectrum. 
This decay was calculated in [52], and again for com¬ 
pleteness we include it here as it was presented in pSEj : 


dN^ 1 _ Qem 1 
dx -1 3-7T X-i 

where r = TOg/m^ and 



) In - + 
r 



T_i(x) =(1 - x)(3 - 2x + 4x 2 - 2x 3 ) 

TT . . „ I? 23 101 2 55 3 „ 

U- 1 {x)=(l-x)\-— + —x-—x+—x 3 (A4) 

+ (3 — 2x + 4x 2 — 2x 3 ) ln(l — x)) 


Note the subscript —1 here is used to remind us this is 
the spectrum calculated in the rest frame of the muon. 
To then obtain the 0-step cascade we would have to ap¬ 
ply Eq. a once, assuming = 2m^/mi -C 1 , and then 
combine this with the FSR spectrum in Eq. |A2| 


3. Annihilations to t + t 

For the case of final state taus, FSR will now be a sub¬ 
dominant contribution. Instead the spectrum will have a 
much larger contribution from leptonic and semi-leptonic 
tau decays: r~ —> v T l~ vi and v T du. The quarks will then 


3. Despite the simplifications afforded, results in this 
framework can be used to estimate the results even 


for general e», as described in Sec. IV 


In this appendix we show how the kinematics of scalar 
cascade decays lead to an expression for the n-step spec¬ 
trum in terms of the (n — l)-step result. In addition we 
outline how Eq. [5] emerges in the small e limit, with error 
0(e 2 ), as well as how the transition to the degenerate 
case as e —> 1 occurs. 

Our starting point is the 0-step spectrum dN^/dx 0 
where xq = 2Eq/itii and Eq is the photon energy in the 
rest frame of </>i. This results from the process <pi -A 'yX, 
where the identity of X depends on the final state con¬ 
sidered. From here we want to calculate dN^/dx 1 - the 
spectrum from a cascade that includes <f >2 —> 4>i(f>i and 
so is one step longer - where x\ = 2 Ei/rri 2 and E\ is 
the photon energy in the <f >2 rest frame. If we assume 
isotropic scalar decays, then we can obtain this by sim¬ 
ply integrating the 0 -step result over all allowed energies 
and emission angles: 


diV 7 
dx 1 


=2 



d cos 6 



dx 0 


dV 7 

dx 0 


5 


^ 2X\ — Xq — COS 9xq 



(Bl) 


where 6 is defined as the angle between the photon mo¬ 
mentum and the <f> 1 boost axis as it is measured in the <j> 1 
rest frame. The limits of integration 0 < Xq < 1 reflect 
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the fact that the photon energy in the 0 i rest frame can 
be arbitrarily soft on the one side, and on the other it 
can have an energy at most half the mass of the initial 
particle, mi/2 here. The S function is simply enforcing 
how the photon energy changes when we move from the 
0i to the 0 2 rest frame, i.e. from Eq to E\. It also sets 
the kinematic range for aq, which is: 


0 < Xi < — ( 1 + 


1 — e i 


(B2) 


Now if we then use the S function to perform the angular 
integral, the one step spectrum reduces to: 


dN 7 _ r 1 '” 1 

dx i ~ 

where we have introduced: 


dx o 


dAA 


Xq^/ 1 — ef dxo 


(B3) 


h,max = min 


2xi 

/l 


I + A/I-Cf 






(B4) 


The maximum here is either set by the maximum physical 
value of Xq, which is 1 , or alternatively by where the <5 
function loses support. We can then repeat this process 
to recursively obtain the ?’th order spectrum from the 
(* — l)th order result. Explicitly we find: 


= 2 r- m 

dxi J t . . 

where we have defined: 


dxi-i 


dNry 


Ci-\\J 1 — dxi —1 


(B5) 


max — min 


2Xi 


i — 1 


2 i_ 


fc=l 


En (i+v 1 -- 


4 1. 




2 Xi 

7T 


l- \/l-e: 


(B 6 ) 


and now the kinematic range of Xi is 


o<* i <-n(i+ 


k =1 


i-4 


(B7) 


With the exact result of Eq. B5 we can now see that in 
the small e limit the result reduces to Eq. [5] with correc¬ 
tions at most of order e 2 , as claimed. The exact result 
also captures an additional feature that the large hier¬ 
archies result does not: the emergence of a degenerate 
step in the cascade as ej —» 1 for some i. As discussed 
in Sec. |IV[ when this occurs, just from the kinematics 
we can see that the (i + l)-step result will reduce to the 
i-step spectrum, b ut sh ifted in energy and normalisation. 

setting 1 — ef = z and then tak- 


B5 


Starting with Eq. 
ing z —} 0 it is straightforward to confirm that the exact 
result also reproduces this behaviour. 


As discussed in Sec. IV there should be a smooth in¬ 
terpolation between the two extreme cases of Cj = 0 and 
£j = 1, and using Eq. |B5| we can demonstrate that indeed 
there is. This is shown in Fig. [lOj where we take the case 
of a 1-step cascade for final state taus with e T = 0.1. We 
plot the two extreme cases and show how intermediate 
e transition between these by plotting five values: 0.3, 
0.5, 0.7, 0.9 and 0.99. Note that as claimed earlier, the 
transition is roughly quadratic in e; for small and inter¬ 
mediate values of e, the result is well approximated by 
the e = 0 result, again highlighting the utility of the large 
hierarchies approximation. 


Appendix C: Model-Building Considerations 


1. A Simple Model 


Let us extend the usual Higgs Portal [13 IMS model to 
include a rich dark sector with n scalar mediators and a 
set of n X 2 symmetries^] This will serve as an illustrative 
example of how different observable signatures depend on 
different model parameters, as discussed in the main text. 

Consider the potential: 


V( X ,<t>i,H) = V x +V H + c k <t>l\H\ 


E 

i=l 


^4 ,i 


2 X <>. - 2"> 




+ 


E 

*)j= 1 


xvy 

4! 


: 


(Cl) 


Here V x and Vh contain the usual mass and quartic terms 
for the DM and Higgs fields. As discussed previously it is 
reasonable that the dark sector is secluded such that the 
dominant portal coupling is Ck4>\\H\ 2 . Upon electroweak 
and Z 2 symmetry breaking the couplings allow an¬ 
nihilations XX *■ 'fii't’i- We assume that DM annihilates 
preferentially to the heaviest mediator through A 4 jn x 2 <(> 2 . 
So it is A 4 that dominantly controls the thermal annihi¬ 
lation cross-section and therefore the DM relic abundance 
U x ft 2 ~ 0.11. The dark sector quartic term will generate 
interactions of the form A ij ((f>i) (j>i(f>j, allowing the media¬ 
tors to cascade decay in the dark sector. Additionally the 
Higgs Portal interaction will generate a mixing between 
and the Higgs. The end result will be a dark cas¬ 
cade ending in the Cfc suppressed decay 0 i —> //, with a 
subsequent photon spectrum that can be fit to the GCE. 

While the thermal relic cross-section depends on A 4 jn , 
the direct detection cross-section will also depend on the 
portal coupling Cfc. This additional small parameter gives 
us the needed freedom to explain the GCE while allevi¬ 
ating constraints from direct detection. Additionally we 
point out that the size of the couplings Ay will need to 


A more complex symmetry structure could allow off-diagonal 
couplings between the scalars and the Higgs, with potentially 
rich observational signatures. We thank Jessie Shelton for this 
observation. 
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be large enough such that decays of the new light states 
occur before BBN. Given the number of new free param¬ 
eters, this setup should not be difficult to construct. Fi¬ 
nally we point out that the Higgs Portal interaction also 
contains a coupling which leads to the decay h —> <f>i<pi. 
Invisible Higgs decay is constrained by collider searches 
which impose an upper bound of about c*, < 10 - 1 2 * 4 5 6 7 8 9 10 11 12 13 m\- 


2. The Sommerfeld Enhancement 

We have seen that the preferred cross-section steadily 
increases with the number of steps in the cascade, moving 
away from the thermal relic value that is favored for the 
direct case. This increased cross-section is also accom¬ 
panied by an increase in the preferred mass scale for the 
DM (indeed, the requirement for a larger cross-section is 
largely driven by the reduced number density of heavier 
DM). In the presence of a mediator much lighter than 
the DM, exchange of such a mediator could enhance the 
present-day annihilation cross-section via the Sommer¬ 
feld enhancement (e.g. [55U59] ). naturally leading to an 
apparently larger-than-thermal annihilation signal. 

However, there are some obstacles to such an inter¬ 
pretation, at least in the simple case we have stud¬ 
ied where the particles involved in the cascade are all 
scalars. For the case of fermionic DM coupled to a light 
scalar or vector of mass with coupling an, the Som¬ 
merfeld enhancement at low velocity is parametrically 
given by m^/anm x . A large enhancement thus requires 
a D > m^/m x . In order to obtain the correct relic den¬ 
sity, we typically require an to be 0(0.01), and so a sig¬ 
nificant Sommerfeld enhancement would require the first 
step in the cascade to involve a mass gap of two orders of 


magnitude. This may be plausible for the electron and 
even muon channels, but is challenging for final states 
involving heavier particles such as taus and 6-quarks; if 
the mediator is heavy enough to decay to these particles, 
the required DM mass becomes much too large to fit the 
GCE even for a one-step cascade, and adding more hier¬ 
archical steps only exacerbates the self-consistency issue 
(as discussed in Secs. m - 


Furthermore, if the DM is a fermion, its annihilation 
into scalars is generically p-wave suppressed, making it 
difficult to obtain a large enough cross-section to ob¬ 
tain the GCE. If instead the DM is a heavy (singlet) 
scalar, the simplest way to couple it to the light scalar 
to which it annihilates is an interaction of the form 
£quartic = When the light scalar obtains a vac¬ 

uum expectation value, this gives rise to an interaction 
of the form A 4 (</>„) </>„x 2 , and repeated exchanges of the 
light scalar tp n can give rise to enhanced annihilation. 
However, assuming (<j> n ) ~ m n , the size of the coupling 
is suppressed by the small mass of the light scalar, even as 
its range is enhanced. Accordingly, a large enhancement 
to annihilation is not expected, at least in this simple 
scenario. 


As discussed in Sec.|IV[ our results can be extended to 
cascades including particles other than scalars, in which 
these later issues do not arise; for example, in the axion 
portal [33], two-step cascades occur through XX sa > 
s —^ aa, a —^ f /*, where s is a dark scalar and a a dark 
pseudoscalar. This annihilation channel is s-wave and 
can be Sommerfeld-enhanced by exchange of the s. How¬ 
ever, the first difficulty described above may still apply, 
with the large hierarchy between the \ and s potentially 
implying a DM mass too large to easily fit the GCE. 
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